home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / os2 / adaptor.zip / ADAPT.ZIP / adaptor / examples / f77 / fftold.f < prev    next >
Text File  |  1993-06-28  |  4KB  |  181 lines
  1.       Program Fast_Fourier_Transform
  2.  
  3. C     Programm zur Berechnung der Potenzen der Einheitswurzeln
  4.  
  5. C     N  Problemgroesse; 2er Potenz
  6.  
  7. C     N2 Anzahl der zu berechnenden Potenzen; N2 = N/2
  8. C     Direkt Fouriersynthese bzw. Fourieranalyse (1/-1)
  9. C     W(N) Feld zur Aufnahme der Potenzen der Einheitswurzeln
  10.  
  11.  
  12.       integer N,N2,k
  13.  
  14.       complex W(:)
  15. cmf$  layout W(:serial)
  16.       complex C(:)
  17.       real F(:), maxf
  18.  
  19.       real Direkt
  20.  
  21.       print * ,'Input k : (N = 2**k)'
  22.       read * ,k
  23.       N = 2**k
  24.       N2 = 2**(k-1)
  25.      
  26.       allocate(W(0:N2-1))
  27.       allocate(C(0:N-1))
  28.       allocate(F(0:N-1))
  29. c
  30. c     random initialization
  31. c
  32.       call cmf_random (f)
  33. c
  34. !HPF$ INDEPENDENT, LOCAL_ACCESS
  35.       do i=0,N-1
  36.         c(i) = f(i)
  37.       end do
  38. c
  39.       Direkt = -1
  40.       call Presort (N,k,C)
  41.       call Roots (W,N2,N,Direkt)
  42.       call Transform(C,W,N,N2,k)
  43. c
  44. c     inverse fft
  45.       Direkt = 1
  46.       call Presort (N,k,C)
  47.       call Roots (W,N2,N,Direkt)
  48.       call Transform(C,W,N,N2,k)
  49. !HPF$ INDEPENDENT, LOCAL_ACCESS
  50.       do i=0,n-1
  51.         f(i) = abs (c(i)/n - f(i))
  52.       end do
  53.       maxf = 0.0
  54. !HPF$ INDEPENDENT, LOCAL_ACCESS
  55.       do i = 0, N-1
  56.          reduce (maxval, maxf, f(i))
  57.       end do
  58.       print *,'Max = ', maxf
  59.       deallocate(F)
  60.       deallocate(C)
  61.       deallocate(W)
  62.  
  63.       end
  64.  
  65. C     *******************************************************
  66.  
  67. C     -------------------------------------------------------
  68. C     ------------ Umordnung der Speicherplaetze ------------
  69. C     ------------     nach Cooley und Tukey     ------------
  70. C     -------------------------------------------------------
  71.  
  72.       subroutine Presort (N,k,C)
  73.       complex C(0:N-1)
  74.       complex H(0:N-1)
  75.       Integer kl,kn
  76.       integer hi, hj, i, index(0:n-1), ix
  77.  
  78. c     initialization of index
  79.  
  80. !HPF$ INDEPENDENT, LOCAL_ACCESS
  81.       do i = 0, n-1
  82.          index(i) = i
  83.       end do
  84.  
  85. c     executing bit reverse for all elements of index
  86.  
  87.       do L=1,k
  88.         kl = 2**(L-1)       
  89.         kn = N / (kl + kl)    
  90.  
  91. !HPF$ INDEPENDENT, LOCAL_ACCESS
  92.         do i = 0, n-1
  93.            ix = index(i)
  94. c          hi = ix / 2
  95.            hi = ishft (ix,-1)
  96. c          hj = hi / kn
  97.            hj = ishft (hi,L-k)
  98. c          hi = mod (hi,kn)
  99.            hi = iand (hi,kn-1)
  100.            index (i) = hi + 2*kn*hj + mod(ix,2)*kn
  101.         end do
  102.        end do
  103.        H = C(index)
  104. !HPF$ INDEPENDENT, LOCAL_ACCESS
  105.        do i = 0, n-1
  106.           C(i) = H(i)
  107.        end do
  108.        end 
  109.  
  110. C     ********************************************************
  111.  
  112.  
  113. C     -------------------------------------------------------
  114. C     ----- Berechnung der Potenzen der Einheitswurzeln -----
  115. C     ----- der Alg. spiegelt eine Rekurrenz erster     -----
  116. C     ----- Ordnung wieder.                             -----
  117. C     -------------------------------------------------------
  118.  
  119.       subroutine Roots (W,N2,N,Direkt)
  120.  
  121.       complex W(0:N2-1),omega
  122. cmf$  layout W(:serial)
  123.       real phi
  124.       real Direkt
  125.  
  126.       phi = Direkt * 8.0 * Atan(1.0) / real(n)
  127.       omega = cmplx(cos(phi),sin(phi))
  128.      
  129.       W(0) = 1
  130.       do i = 1,N2-1
  131.          W(i) = W(i-1) * omega
  132.       end do
  133.       end 
  134.  
  135.  
  136.  
  137. C     ********************************************************
  138.  
  139.  
  140. C    **********************************************************
  141.  
  142.  
  143. C    --------------------------------------------
  144. C    ----- Transformation nach Cooley-Tukey -----
  145. C    --------------------------------------------
  146.  
  147.       subroutine Transform(C,W,N,N2,k)
  148.  
  149. C     Feld c repraesentiert die Fourierkoeffizienten
  150. C     Feld w beinhaltet die Potenzen der Einheitswurzeln
  151.  
  152.       complex C(0:N-1),W(0:N2-1),H(0:N-1)
  153. cmf$  layout W(:serial)
  154.       complex t
  155.  
  156.       integer kl,kn,kl1,hi,index(0:N-1),ni
  157.       
  158.       do L=1,k
  159.         kl = 2**(L-1)
  160.         kn = N2/kl
  161.         kl1 = 2 * kl
  162. !HPF$ INDEPENDENT, LOCAL_ACCESS
  163.         do i = 0, n-1
  164.            index(i) = ieor (i,kl)
  165.         end do
  166.         H = C(index)
  167. !HPF$ INDEPENDENT, LOCAL_ACCESS
  168.         do i = 0, n-1
  169. c          hi = mod (i,kl1)
  170.            hi = iand (i,kl1-1)
  171. c          ni = mod (i,kl)
  172.            ni = iand (i,kl-1)
  173.            if (hi .ge. kl) then
  174.               C(i) = H(i) - C(i) * w(ni*kn)
  175.             else
  176.               C(i) = C(i) + w(ni*kn) * H(i)
  177.            end if
  178.         end do
  179.       end do
  180.       end 
  181.